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1  General  information 

This  is  the  final  technical  report  on  research  carried  out  under  Air  Force  Grant  No. 
F49620-92-J-0158-DEF,  in  the  period  April  15,  1992  -  April  14,  1995.  The  project  ti¬ 
tle  was  “Local  preconditioning  of  the  Euler  equations  and  its  numerical  applications.” 
It  was  awarded  for  research  exploiting  the  recent  discovery  of  local  preconditioning 
matrices  that  minimize  the  spread  among  the  wave  speeds  admitted  by  the  Euler 
equations  for  any  Mach  number.  Closely  related  to  this  research  is  the  work  carried 
out  under  Augmentation  Grant  No.  F49620-93-1-0417;  its  project  title  is  “Efficient  ex¬ 
plicit  integration  schemes  for  the  hyperbolized  Navier-Stokes  equations.”  This  project 
considers  the  application  of  local  preconditioning  to  a  stiff  hyperbolic  system  describ¬ 
ing  viscous  conducting  flow  with  a  finite  relaxation  time.  Over  the  past  year  part  of 
the  effort  in  this  subject  was  diverted  to  the  main  subject  in  order  to  help  solve  the 
problem  of  the  stagnation-point  instability  (see  below),  which  was  halting  progress 
in  l)otli  principal  and  augmentation  projects.  Where  relevant,  research  findings  from 
the  augmentation  project  are  included  below. 

2  Research  in  the  past  year 

Since  the  previous  report  research  under  the  principal  grant  has  focussed  on  the 
problem  of  instability  in  stagnation  regions,  experienced  in  the  use  of  a  variety  of 
inviscid  flow  codes:  conservative  and  nonconservative,  cell-average  and  cell-vertex 
based. 

In  July-.'Vugust  '95  a  concentrated  effort  (not  funded  by  AFOSR)  of  some  develop¬ 
ers  and  users  of  local  preconditioning  at  ICASE/NASA  Langley  Research  Center  led 
to  a  modification  of  the  Van  Leer  preconditioning  that  makes  this  matrix  insensitive 
to  the  flow  direction  near  stagnation  points,  where  the  flow  direction  is  ill-determined. 
The  modified  matrix  immediately  cured  the  stability  and  accuracy  problems  for  the 
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cell-vertex,  fluctuation- split,  nonconservative  code  of  Lisa  Mesaros  (doctoral  candi¬ 
date,  U.  of  MI);  the  conservative  version  of  her  code  was  not  helped  very  much  by  the 
modification.  This  suggested  that  there  was  yet  another  cause  of  the  stagnation-point 
difficulties. 


Subsequently,  David  Darmofal  (NSF/AFOSR  post-doc,  U.,  of  MI),  in  collabora¬ 
tion  with  Peter  Schmid  (U.  of  WA),  studied  the  eigenvector  structure  of  the  matrices 
arising  in  the  preconditioned  Euler  equations.  He  concluded  that  the  main  problem 
for  vanishing  flow  speed  or  Mach  number  M  is  the  large  deviation  from  orthogonality 
of  the  eigenvectors:  certain  eigenvectors  become  parallel  for  M  -y  0.  When  marching 
to  a  steady  solution  this  causes  a  strong  initial  transient  that  may  lead  to  instability. 

Based  on  this  analysis  a  modification  waa  proposed  for  preconditioning  matrices 
of  the  Van  Leer-Turkel  family;  the  modified  matrices  were  successfully  implemented 
in  an  unstructured  conservative  code.  This  marks  the  first  time  that  the  stability 
problem  has  been  removed  in  the  use  of  a  conservative  code.  The  modification  consists 
of  downward  limiting  of  the  Mach-number  value  used  in  the  (1,1)  matrix  element 
(~  M^),  e.g.  by  a  minimum  value  that  is  a  small  fraction  of  the  far-field  Mach  number. 
In  previous  “fixes”  the  Mach  number  was  limited  everywhere  in  the  preconditioning 
matrix  and  its  products  with  other  matrices;  apparently  this  does  not  affect  the 
numerical  fluxes  in  a  desirable  way. 

We  expect  the  new  modified  preconditioner  to  solve  the  stagnation  problems  of 
other  conservative  preconditioned  codes;  this  is  currently  under  investigation  in  other 
projects. 

The  application  of  the  Van  Leer  preconditioning  to  the  design  of  multi-stage 
marching  scheme  with  optimal  high-frequency  damping,  and  the  use  of  such  schemes 
in  a  multi-grid  relaxation  framework,  are  the  subject  of  a  thesis  by  John  Lynn,  de¬ 
fended  in  our  Department  on  May  17,  1995.  He  has  shown  that  the  double  benefit 
of  preconditioning,  known  from  one-dimensional  calculations  by  Tai  (Ph.D.,  U.  of 
MI,  1990),  persists  in  two  dimensions.  The  preconditioning  already  speeds  up  single¬ 
grid  relaxation;  the  guaranteed  high-frequency  damping  of  the  optimized  multi-stage 
schemes  provides  an  extra  speed-up  of  a  factor  2-3,  and  additional  robustness. 

The  thesis  includes  tables  of  recommended  multi-stage  coefficients,  and  also  dis¬ 
cusses  the  extensions  to  Navier-Stokes  operators  and  to  unstructured  grids.  In  the 
numerical  examples  Lynn  avoided  the  stagnation-point  instability  by  avoiding  flow 


problems  including  stagnation  points. 

The  above  research  has  been  reported  at  the  12th  AIAA  Conference  on  Compu¬ 
tational  Fluid  Dynamics,  June  19  -  22,  San  Diego,  CA.  Two  conference  papers  were 
prepared  under  the  principal  AFOSR  grant;  these  are  listed  below. 


1.  B.  van  Leer,  L.  Mesaros,  C.-H.  Tai  and  E.  Turkel,  “Local  preconditioning  in  a 
stagnation  point,”  AIAA  Paper  95-1654CP,  1995. 
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2.  J.  F.  Lynn  and  B.  van  Leer,  “A  semi-coarsened  multi-grid  solver  for  the  Eu¬ 
ler  and  Navier- Stokes  equations  with  local  precnditioning,”  AIAA  Paper  95- 
1667CP,  1995. 

The  remainder  of  the  research  was  dedicated  to  further  development  of  the  pre¬ 
conditioning  matrix  for  the  Navier-Stokes  equations  so  as  to  make  it  useful  for  very 
low  cell  Reynolds-numbers.  Preconditioning  in  the  presence  of  a  one-equation  turbu¬ 
lence  model  (Spalart-Allmaras)  has  also  been  investigated.  Neither  subject  has  been 
rounded  off. 


3  Review  of  results  from  entire  contract  period 

During  the  contract  period  of  three  years  major  progress  was  made  in  the  under¬ 
standing,  application  and  improvement  of  local  preconditioning  for  the  Euler  and 
Navier-Stokes  equations.  One  of  the  foremost  results  is  the  demonstration  that  the 
use  of  our  Euler  preconditioner  in  two  space  dimension  yields  a  double  benefit  in  a 
multi-grid  strategy;  a  single-grid  speed-up  factor,  owing  to  removal  of  system  stiffness, 
and  an  independent  multi-grid  speed-up  factor,  owing  to  guaranteed  high-frequency 
damping.  These  results  are  contained  in  the  Ph.D.  thesis  of  John  Lynn,  whose  doc¬ 
toral  work  was  supported  by  the  contract  over  the  last  18  months.  It  is  worth  noticing 
that  among  researchers  that  use  local  preconditioning  for  the  sake  of  multi-grid  re¬ 
laxation,  almost  all  select  point-Jacobi  preconditioning,  which  only  gives  the  second 
bewnefit. 

Basic  understanding  of  the  potential  possibilities  present  in  the  family  of  optimal 
Euler  preconditioners  was  obtained.  It  turned  out  that,  in  spite  of  the  overwhelming 
number  of  free  parameters,  much  of  the  parameter  space  is  irrelevant.  The  symmetric 
Van  Leer-Lee-Roe  2  parameter  preconditioning  matrix,  on  which  this  contract  was 
based,  turns  out  to  be  the  most  desirable  form  under  all  circumstances.  For  very 
low  Mach  number,  it  offers  the  flexibility  to  have  the  condition  number  increase  by 
a  factor  2  in  order  to  reduce  flow-angle  sensitivity.  The  matrix  has  also  caused  a 
break-through  in  the  formulation  of  genuinely  multi-dimensional  fluctuation-splitting 
schemes,  after  it  was  discovered  that  it  yields  the  proper  decoupling  of  the  acoustic 
from  the  convective  components  present  in  the  Euler  equations. 

Robustness  of  Euler  calculations  in  the  presence  of  preconditioning  has  been  an 
issue  during  the  entire  contract  period,  but  it  appears  that  the  worst  problems  are 
over,  now  that  the  stagnation-point  instability  has  been  properly  understood.  This 
instability  is  caused  in  part  by  the  lack  of  orthogonality  of  certain  matrix  eigenvectors, 
resulting  in  a  transient  up-swinging  of  numerical  errors. 

The  Euler  preconditioning  has  been  embedded  in  a  Navier-Stokes  preconditioning 
through  the  use  of  a  simple  Jacobi  block  to  account  for  the  dissipative  terms.  There  is 
still  a  robustness  issue  in  the  limit  of  vanishing  cell  Reynolds-number.  The  inclusion 
of  a  convective  turbulence-transport  equation  appears  to  be  easily  accounted  for  in  the 
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preconditioning  matrix,  in  the  same  way  as  is  done  for  the  extra  continuity  equations 
needed  in  describing  multi-species  reacting  flow. 

All  results  not  included  in  Lynn’s  thesis  will  be  presented  in  the  doctoral  thesis 
of  Dohyung  Lee,  which  is  expected  to  be  defended  in  the  Winter  Term  of  1996. 

When  reviewing  progress  in  CFD  over  the  past  three  years,  it  is  observed  that 
the  introduction  of  the  Van  Leer-Lee-Roe  preconditioning  matrix  at  the  10th  AlAA 
CFD  Conference  in  June  1991  set  a  new  CFD  trend.  Local  preconditioning  matrices 
for  the  low-speed  flow  regime  had  been  introduced  by  Turkel  from  1984  onwards, 
and  were  used  extensively  by  Merkle  et  al.  at  Pennsylvania  State  University,  but 
their  use  did  not  spread  further.  At  present  there  is  a  rapidly  growing  appreciation  of 
local  preconditioning,  both  in  the  USA  and  Europe  (Netherlands,  Belgium,  Germany, 
UK,  France),  probably  because  the  locality  of  this  convergence-acceleration  technique 
combines  well  with  parallel  computer  architectures. 

4  Future  research 

To  make  our  type  of  local  preconditioning  a  robust  technique,  as  robust  as  the  tradi¬ 
tional  point-Jacobi  preconditioning,  and  suited  for  non-expert  users,  some  work  still 
needs  to  be  done: 

a.  Testing  the  robustness  of  the  stagnation-point  modifications  under  a  wide  range 
of  circumstances,  i.  e.  Mach  numbers,  flow  angles,  cell  aspect-ratios  and  cell 
Reynolds-numbers; 

b.  Optimizing  the  smoothing  of  matrix  elements  in  the  sonic-point  singulatity; 

c.  Overcoming  the  loss  of  robustness  for  very  low  cell  Reynolds  numbers. 

Among  new  uses  of  local  preconditioning  awaiting  exploration  are: 

•  Local  preconditioning  for  atmospheric  problems,  eliminating  stiffness  due  to 
acoustic  and  gravitational  waves; 

•  Preconditioning  for  magneto-hydrodynamics,  eliminating  stiffness  due  to  magneto- 
acousic  and  Alfven  waves; 

•  Local  preconditioning  expanded  hydrodynamic  equation  systems  derived  from 
the  collisional  Boltzmann  equation  (a  continuation  of  the  work  under  the  current 
augmentation  grant).  This  has  applications  in  upper-atmosphere  and  plasma 
flows  as  well  as  micro-manufacturing  and  solid-state  modeling. 

We  plan  to  submit  a  new  research  proposal  addressing  these  topics  in  the  near  future. 
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Abstract 

Local  preconditioning  yields  desirable  effects 
such  as  convergence  acceleration  and  preser¬ 
vation  of  accuracy  in  the  incompressible  limit, 
but  these  come  at  the  expense  of  robust- 
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ness.  Current  research  focuses  on  making  the 
various  preconditioned  Euler  discretizations 
that  have  been  developed  more  reliable  near 
flow  singularities  such  as  sonic  and  stagnation 
points.  In  this  report  we  discuss  the  loss  of 
stability  in  computing  stagnating  flow  with 
the  symmetric  Van  Leer-Lee-Roe  precondi¬ 
tioning,  and  do  an  exhaustive  search  for  modi¬ 
fications  of  the  matrix  (for  low  Mach  numbers) 
by  which  the  instability  can  be  avoided.  Some 
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numerical  support  regarding  the  effectiveness 
of  these  modifications  is  presented. 


1  Introduction 


Local  preconditioning  can  be  used  to  make  the 
time-dependent  Euler  equations  more  suited 
to  the  numerical  computation  of  steady  solu¬ 
tions.  In  the  following  analysis  precondition¬ 
ing  is  accomplished  by  multiplying  the  spa¬ 
tial  differential  operator  by  a  locally  evaluated 
matrix.  First,  let  us  introduce  some  notation. 

The  two-dimensional  Euler  equations  will 
be  written  as 


dU 

dt 


=  0, 


(1) 


where  the  state  variables  axe  defined  differen¬ 
tially: 

dU  =  {—,du,dv,dSf-,  (2) 

pa 

a  denotes  sound  speed,  S  entropy.  These  vari¬ 
ables  have  the  property  of  symmetrizing  the 
Euler  equations: 
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(4) 


Much  of  the  analysis  can  be  done  assum¬ 
ing  that  the  fluid  moves  in  the  positive  x- 
direction,  so  that  u  equals  the  full  flow  speed  q, 


and  V  vanishes  (but  not  its  derivatives).  The 
Euler  equations  then  read 


dU 

dt 
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dU 

=  0,  (5) 
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0 

0 

0 

dy 

0 

0 

0 

where  M  is  the  Mach  number. 

The  preconditioned  Euler  equations  read 

where  P{U)  is  a  positive  definite  matrix.  If  P 
is  symmetric,  so  is  and  the  system 

P(t/)-f +  +  =  0  (7) 

is  symmetric,  just  as  the  original  system.  If 
P  is  not  symmetric,  it  is  desirable  that  the 
preconditioned  system  be  symmetrizable;  th 
puts  certain  constraints  on  the  elements  of  P 
[11- 

In  the  course  of  our  research  on  local  pre¬ 
conditioning  of  the  Euler  equations  the  follow¬ 
ing  computational  benefits  have  been  identi¬ 
fied. 


1.  Local  preconditioning  can  be  designed  to 
remove  the  stiffness  of  the  system  of  equa¬ 
tions  caused  by  the  range  of  the  charac¬ 
teristic  speeds,  thus  improving  the  con¬ 
vergence  rate  of  any  discrete  marching 
scheme  [2]. 
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2.  It  then  makes  the  system  of  equations  be¬ 
have  more  like  a  scalar  equation,  which  is 
advantageous  in  designing  and  applying 
additional  convergence-acceleration  tech¬ 
niques  [3]. 

3.  It  generates  spatial  discretizations  that 
may  retain  their  accuracy  at  low  Mach 
number,  in  contrast  to  standard  dis¬ 
cretizations  [4,  5,  6]. 

4.  It  can  be  designed  to  decouple  the  acous¬ 
tic  equations  from  the  convective  equa¬ 
tions,  allowing  genuinely  multidimen¬ 
sional  discretizations  [1] 


These  desirable  properties  come  at  the  ex¬ 
pense  of  robustness.  Current  research  focuses 
on  making  the  various  preconditioned  Euler 
discretizations  that  have  been  developed  more 
reliable  near  flow  singularities,  such  as  sonic 
and  stagnation  points.  In  this  report  we  de¬ 
scribe  the  loss  of  stability  in  computing  stag¬ 
nating  flow  with  the  symmetric  Vau  Leer- Lee- 
Roe  preconditioning  [2],  and  how  to  overcome 
it  by  modifying  the  matrix  for  low  Mach  num¬ 
bers. 


2  Sensitivity  to  flow  angle 

With  the  choice  of  flow  variables  (2)  and  the 
flow  aligned  with  the  x-axis,  the  Van  Leer- 
Lee- Roe  preconditioning  matrix  for  subsonic 


flow  (M  <  1)  takes  the  form 


m)  = 


I  Ml  _M 


M 
"  P 
0 

0 


0 

0 


0  0  \ 

0  0 

13  0 

0  1/ 


(8) 


where  13  =  \/l  —  M^.  Use  of  this  matrix  re¬ 
duces  the  condition  number  of  the  character¬ 
istic  speeds  from  (M  -f- 1)/  min(l  —  M,  M)  to 
l/y5.  For  an  appreciation  of  the  stability  prob¬ 
lem  near  a  stagnation  point  we  must  obtain 
the  form  of  the  matrix  valid  for  any  flow  an¬ 
gle  (i>,  i.  e. 

pm  =  B:;\u)P{U)Rm  (9) 


with 


Rm  = 


/I  0  0  0  \ 

0  cos  (f>  sin  0 

0  —  sin<^  cos^  0  ’ 

\  0  0  0  1  / 


(10) 


this  yields  Equation  (11).  It  is  seen  that 
this  matrix  remains  sensitive  to  the  flow  an¬ 
gle  when  the  Mach  number  decreases  toward 
zero,  because  its  inner  elements  (2,2),  (2,3), 

(3.2)  and  (3,3),  which  depend  on  <i>,  remain 
0(1).  In  this  case  numerical  perturbations 
to  u  and  v  that  are  small  in  absolute  value 
may  not  be  small  when  compared  to  the  values 
of  u  and  v,  causing  0(1)  variations  in  (j)  and 
the  four  matrix  elements.  This  sensitivity  can 
be  traced  back  to  the  fact  that  the  elements 

(2.2)  and  (3,3)  of  P  are  not  equal.  It  is  be¬ 
lieved  to  contribute  to  the  numerical  instabil¬ 
ity  near  a  stagnation  point,  often  experienced 
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pm  = 


/ 
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_Mcos<^ 

0 


—  ^  cos  <l> 
cos^  ^  sin^  <l> 
sin  <j)  cos  (j) 

0 


—ysm<t>  0  ^ 

+  1  —  sin  (j)co3<t>  0 

^  SlXp  (j)  +  ^  cos^  <l>  0 

0  1  / 


(11) 


with  conservative  upwind  Euler  schemes  that 
include  the  above  matrix.  In  such  schema 
the  artificial-dissipation  matrices  [2]  used  in 
the  streamwise  and  normal  fluxes  are 


P-'IPAI,  P'^\PBl 


(12) 


eral  elements  of  the  matrix  (8)  near,  a  stag¬ 
nation  point,  especially  the  (1,1)  element.  In 
practice  it  has  been  found  necessary  to  bound 
the  value  of  M  away  from  zero  in  these  el 
ments,  e.g.  replace  it  by 


rather  than  1A|  and  \B\,  as  in  the  original 
upwind  schemes.  In  the  update  the  spatial 
residual  in  each  cell  is  multiplied  by  the  cell- 
centered  value  of  P^,  creating  products 


Af  =  min(0.1A!foo?0-l)'  (I'l^) 


(  )  center  (  P<j>  ) 


(13) 


that  may  vary  erratically  near  a  stagnation 
point  and  deviate  appreciably  from  the  v^- 
ues  elsewhere  in  smooth  flow,  which  should  be 
close  to  the  identity  matrix.  In  a  study  by  D. 
Lee  17]  a  uniform  slow  flow  perturbed  by  ro¬ 
tating  the  velocity  in  one  cell  by  90®  becanue 
unstable  when  advanced  in  time  by  explici 
first-order  upwind  scheme  preconditioned  as 
above.  The  instability  could  be  forestalled, 
but  not  avoided,  by  taking  smaller  time  steps. 
A  similar  behavior  was  found  when  simulating 
an  isolated  stagnation  region.  Andrew  God¬ 
frey  reports  [8]  that  implicit  time  integration 
can  suppress  the  instability  if  the  grid  used  is 
not  too  fine. 

Besides  flow-angle  sensitivity  there  is  a 
more  obvious  problem:  the  vanishing  of  sev- 


For  conservative  schemes,  which  include  the 
inverse  of  P  in  the  artificial-viscosity  coeffi¬ 
cients  (12),  such  a  procedure  clearly  makes 
sense  [2];  less  obvious  is  it  that  nonconser¬ 
vative  schemes  also  increase  in  robustness  by 
this  measure.  It  appears  that  use  of  the  ac¬ 
tual  value  of  M  makes  the  pressure  correc¬ 
tion  allowed  by  the  preconditioned  equations 
so  small  that  the  flow  does  not  properly  tur- 

in  the  stagnation  region. 

Returning  to  the  subject  of  flow-angle  sen- 
sitiviy,  an  indication  that  this  indeed  con 
tributes  to  stagnation-region  instability  was 
found  by  Tai  [9],  who  obtained  stability  and 
convergence  by  using  a  nonconservahve  ver¬ 
sion  of  the  preconditioned  schemes,  m  which 
the  product  (13)  is  omitted. 
preconditioned  schemes  based  on  Turkel  s  [lOJ 
matrix,  whether  or  not  in  conservation  form, 
do  not  exhibit  the  instability.  This  can  be 
readily  understood  from  the  structure  of  this 
matrix;  for  streamwise  x-axis  and  sufficien  y 
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small  Mach  number  it  reduces  to 


Pt{U)  = 


/  0  0  0  \ 

-M  1  0  0 

0  0  10’ 

\  0  0  0  1  / 


(15) 


which  is  seen  to  be  lower  triangular.  The  ele¬ 
ments  (2,2)  and  (3,3)  of  this  matrix  are  equal; 
in  consequence,  its  version  for  arbitrary  flow 
angle, 


Pt4>{U)  = 


/  0  0  0  \ 
—Mcos<f)  10  0 
—Msm4>  0  10 
0  0  0  1 


(16) 


is  well  behaved  for  M  0.  It  is  possible 
to  extend  the  matrix  (15)  so  as  to  yield  the 
minimum  possible  spread  of  the  characteristic 
speeds  throughout  the  subsonic  domain,  with¬ 
out  giving  up  its  lower-triangular  structure. 
The  resulting  matrix  is 


/ 

Pt{U)  = 

\ 


Mi 

-M 

0 

0 

0 


0  0  0  \ 
1  0  0 
0^0 
0  0  1/ 


(17) 


While  it  does  not  suffer  from  flow-angle  sen¬ 
sitivity,  the  preconditioning  matrix  (15)  has 
another  problem:  it  lies  at  the  limit  of  ad¬ 
missibility.  Turkel  [10]  found  that  the  pre¬ 
conditioned  system  of  equations  is  no  longer 
symmetrizable,  implying,  among  other  things, 
that  it  no  longer  has  an  entropy  condition  [11] 
distinguishing  between  physical  and  nonphys¬ 
ical  admissible  solutions.  An  arbitrarily  small 
perturbation  of  the  matrix,  though,  can  make 


the  system  symmetrizable  again.  Indeed,  it 
has  been  pointed  out  by  Turkel  that  the  ma¬ 
trix  (15)  itself  does  not  work  in  practice,  but 
the  perturbed  matrix 


Pxm  = 


/  {1  +  e)M^  0  0  0  \ 

-M  10  0 

0  0  10’ 

0  0  0  1  / 


(18) 


with  c  >  0  for  symmetrizability,  does  speed 
up  numerical  convergence. 

It  is  not  a  priori  clear  that,  in  the  limit 
of  M  — »  0,  the  matrix  (15)  is  the  only  opti¬ 
mal  preconditioner  with  the  property  that  its 
(2,2)  and  (3,3)  element  axe  equal.  It  would  be 
preferable  if  an  optimal  matrix  existed  closer 
to  Van  Leer’s,  i.  e.  more  nearly  symmetric.  To 
find  out  about  this,  the  best  one  can  do  is  to 
examine  all  possible  matrices.  This  is  actually 
done  in  the  next  section. 


3  Analysis  in  the  incom¬ 
pressible  limit 

Since  the  stability  problem  with  the  symmet¬ 
ric  preconditioner  occurs  only  for  small  M,  it 
is  advantageous  to  base  the  further  analysis  on 
the  incompressible  Euler  equations,  at  least  as 
a  start.  This  means  the  reference  Mach  num¬ 
ber,  e.  g.  Moo,  is  assumed  to  be  small;  the  flow 
velocity  itself  is  now  regarded  to  be  0(1). 

The  Euler  equations  for  incompressible  flow 
are  made  hyperbolic  through  “artificial  com¬ 
pressibility”  [12],  implemented  by  adding  a 
time  derivative  of  pressure  to  the  elliptic  con- 
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tinuity  equation: 


0,  (19) 

0,  (20) 

0,  (21) 

where  a  is  now  a  constant  artificial  speed  of 
sound,  and  p  is  also  a  constant;  the  entropy 
equation  drops  out.  Using  otherwise  the  same 
variables  and  coordinates  as  before,  the  equa¬ 
tions  can  be  written  as 


raipidly  with  m.  Thus,  6,  c,  /  and  F  can  all 
be  0(1),  but  h  and  c  must  be  equal  or  close 
to  equal,  while  /  and  F  must  be  opposite  or 
close  to  opposite.  For  the  present  purpose  it 
therefore  suffices  to  assume 

6  =  c;  (26) 

F  =  -/i  (27) 

it  is  instructive  to  carry  b  and  c  along  as  sepa¬ 
rate  parameters  2is  long  as  possible,  since  this 
allows  us  to  link  the  Vaji  Leer  and  Turkel  m 
trices. 


+  War  +  Vj,  = 

pa^ 

Ut  +  UUx  -b  VUy  +—  = 

P 

Px 

Vt  -f  UVx  -b  VVy  H - = 


dU 

dt 


+  a 


0  1  0  \ 
1  m  0 
0  0  m  y 


dx 


Optimad  asymmetric  matri¬ 
ces 


/  0  0  l\  QTT 

-b  a  0  0  0  ^=0,  (23) 

V  1  0  9  ) 

where  m  =  uj a  is  the  artificial  Mach  number, 
of  magnitude  0(1).  Next  we  precondition  this 
system  with  the  most  general  matrix  possible: 


P{U)  = 


/  a  D  E\ 
d  b  F  y, 

\e  f  c  ) 


(24) 


note  that  a  is  a  free  paj2uneter  and  has  noth¬ 
ing  to  do  with  the  speed  of  sound.  For  an 
arbitrary  flow  angle  this  matrix  transforms 
into  P,j,(U),  given  in  Equation  (25).  Consider 
first  the  2x2  block  of  elements  (2,2),  (2,3), 
(3,2),  (3,3).  In  order  to  avoid  too  strong  a  de¬ 
pendence  on  the  flow  angle  near  a  stagnation 
point,  i.  e.  for  m  — >  0,  both  b  —  c  and  f  F 
must  either  equal  zero  or  vanish  sufficiently 


Next  we  study  the  characteristic  speeds  of  the 
preconditioned  equations,  i.  e.  the  propaga¬ 
tion  speeds  of  the  plane- wave  solutions  admit¬ 
ted  by  these  equations.  If  we  denote  the  prop 
agation  direction  by  the  angle  9,  the  speeds 
are  the  eigenvalues  of  the  matrix 

P{Acos9  +  BsinO).  (2 

Using  the  coefficient  matrices  from  (23)  and 
the  general  preconditioning  (24),  the  eigen¬ 
value  equation  becomes  (29).  or 

-  -b  K2\^  -  KiX  +  Ko  =  0,  (30) 


with 

K2  =  [d  +  D  +  bm  +  cm]  cos  6 

+  {e  +  E)sm9,  (31) 

K\  =  \[d D  +  bm)cm  —  ab dD 
-  eFm  -  fEm  -  fFni^]  cos^  9 
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+  (— ac  +  e£J)  sin^  0 
+  [(/  +  +  (e  +  E)bm  +  {E  -  fm)d 

+  {e  —  Fm)D]sinOcosd,  (32) 

Kq  =  —m  cos  6  det  P.  (33) 

If  P  were  an  optimal  preconditioner,  all  wave 
speeds  would  turn  out  equal,  i.  e.,  the  acoustic 
waves  would  propagate  at  the  flow  speed  in  all 
directions.  This  suggests  the  following  target 
equation  for  the  A’s: 

{mcosO  —  A)(A^  —  m^)  =  0,  (34) 

or 

—  A^  +  A^m  cos  $  +  Am^  —  m®  cos  ^  =  0.  (35) 

Comparing  the  target  equation  to  the  gen¬ 
eral  equation  leads  to  6  constraints  (two  from 
K21  three  from  Ki  and  one  from  Kq)  on  the 
9  parameters  a,  6,  c,  d,  e,  D,  E.  One  therefore 
would  hope  to  find  a  one-parameter  family  of 
matrices  satisfying  the  additional  constraints 
(26)  and  (27). 

Starting  with  the  sin^-term  in  K21  we  see 
that 

e  -I-  E  =  0;  (36) 


it  does  not  bode  well  for  symmetrizability 
that  both  pairs  {e,E)  and  (/,  F)  must  be 
antisymmetric.  Using  both  (27)  and  (36),  the 
remaining  conditions  reduce  to 

d  +  D  +  hm  =  m(l  —  c)  (37) 

from  K2, 

(d  +  D  +  bm)cm  —  ab-\-  dD 
-|-(e  +  frrif  -  ^  -  -m?,  (38) 

ac  -1-  (39) 

and 

{D  -  d)(e  -h  Em)  =  0  (40) 

from  Ki ,  and 

{ab  —  dD)c  -I-  ap  -I-  be^ 

-{d  +  D)ef  =  m‘^  (41) 

from  Kq. 

Considering  Eq.  (40)  it  appears  we  have 
two  choices: 

e  +  fm  =  0,  (42) 

which  also  means 

F  -I-  Fm  =  0,  (43) 
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or 

d  =  D.  (44) 

Pursuing  first  the  choice  e  +  fm  =  0  we  find 
that  Eqs.  (37,38,39,41),  while  up  to  cubic  in 
the  unknown  parameters,  allow  elimination  of 
a,  6,  d  and  D.  The  result  is  a  relation  between 
c  and  /, 

(c2  +  _  I)(c2  +  _  c)  =  0,  (45) 


Going  back  to  Eq.  (40),  we  still  have  to  fol¬ 
low  its  second  branch  Eq.  (44),  which  makes 
the  pair  (d,  £>)  symmetric.  Again  the  alge¬ 
bra  becomes  complicated,  forcing  us  to  insert 
b  =  c  right  away.  When  eliminating  a  and  d 
from  the  remaining  equations  it  turns  out  that 
c  also  drops  out,  yielding 

(e  + /m)^  = -^m^,  (50) 


meaning  either 

f  =  l-c^  (46) 

or 

f  =  c(l  -  c).  (47) 

The  first  choice,  Eq.  (46),  suggests  interpret¬ 
ing  c  and  /  as  cos^’  and  sinV’,  where  rp  is 
an  arbitrary  angle;  it  leads  to  the  following 
two-parameter  matrix  (or  its  transposed)  sat¬ 
isfying  all  previous  constraints  except  (26): 


which  can  not  be  satisfied  by  real  e  and  /. 
Apparently,  the  equality  of  d  and  D  exclude 
the  equality  of  b  and  c. 

Upon  inspecting  the  yield  of  the  above  anal¬ 
ysis,  viz.  matrices  (48)  and  (49),  it  appears 
that  equality  of  b  and  c,  with  c  always  less  than 
1,  makes  these  matrices  even  more  skewed 
than  Turkel’s  matrix,  and  therefore  further 
away  from  producing  a  symmetrizable  precon¬ 
ditioned  system.  The  extra  freedom  offered  by 
the  parameters  e  and  /,  hitherto  not  included 


P(U)  = 


'I  —  c^m 


(1  —  b)m 
b 

±y/r^ 


T  —  c^m 


in  any  analysis,  does  not  lead  to  more  prefer¬ 
able  matrices,  and  we  shall  henceforth  assume 
these  to  be  zero.  For  future  reference  we  re- 
•peat  the  matrix  (48),  inserting  c  =  1,  whii 
eliminates  e  and  /: 


This  general  formula  includes  the  Turkel  ma¬ 
trix  (6  =  c  =  1)  and  the  Van  Leer  matrix 
(6  =  2,  c  =  1);  we  shall  discuss  its  properties 
further  below. 

The  second  choice,  Eq.  (47),  complicates 
the  algebra;  to  simplify  the  analysis  we  have 
to  implement  the  constraint  (26)  right  away. 
With  6  =  c  we  obtain  the  one-parameter  ma¬ 
trix  (49).  In  this  matrix  we  may  replace  c 
by  |(1  -t-cos  <p),  whereupon  /  becomes  ^sintp. 
Turkel’s  matrix  again  is  obtained  for  c  =  1. 


/  (1  —  6)m  0  ^ 

P{U)  =  -m  6  0  ;  (51) 

Vo  0  1/ 

this  is  the  incompressible  version  of  the  fol¬ 
lowing  matrix  valid  for  subsonic  compressible 
flow: 

f  (l-6)f  0  0 
-f  b  0  0 
0  0^0 
0  0  0  1 
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P{U)  = 


(49) 


One  choice  of  6(M)  that  links  the  Turkel  ma¬ 
trix  for  M  =  0  to  the  Van  Leer  matrix  for 
greater  subsonic  values  of  M  is 

b{M)  =  i  +  1  -  (53) 


The  family  of  matrices  of  the  form 


P(U)  = 


a  d  0  0  y 

d  b  0  0 

0  0  c  0  ’ 

0  0  0  1  / 


(56) 


where  r  is  some  positive  power,  e.  g.  r  =  2. 

3.2  Sub-optimal  symmetric  ma¬ 
trices 

We  are  now  left  with  two  possible  strategies 
of  fighting  the  stagnation-point  instability: 

1.  use  a  matrix  of  the  form  (52),  with  b(M) 
varying  with  the  Mach  number  from 


5(0)  =  1 

(54) 

6(1)  =  2. 

(55) 

2.  develop  a  one-paxcimeter  family  of  sub- 
optimal  symmetric  matrices  satisfying 
the  same  conditions  (54)  and  (55). 

In  this  section  we  shall  pursue  the  second  pos¬ 
sibility.  Preserving  symmetry,  even  if  it  means 
giving  up  optimality,  is  valuable  because  it  is 
cruicial  in  coding  update  schemes  based  on 
multidimensional  fluctuation  splitting  [13];  see 
further  Section  4  on  numerical  results. 


has  been  analyzed  in  great  detail  by  Wen- 
Tzong  Lee[4].  From  this  analysis  we  conclude 
that  the  following  variation  on  the  Van  Leer 
matrix  (8)  will  offer  suflScient  freedom: 

/  Ci^  -Cif  0  0  \ 

P(U)=  C'i^  +  C'2  0  0 

^  ^  0  0  040  0  ' 

^  0  0  0  C3  / 

(57) 

With  this  choice  of  preconditioning,  the  plane- 
wave  speeds  A(^)  follow  from  the  equation 

{C2M cose  -  xxCsM cos e  -  X)  X 
{A^ -h  AM(C4 -C'2)cos0  (58) 

— C'iC'4M^[(1  —  M^)  cos^  e  -t-  sin^  ^.]}  =  0; 

for  0  =  0  we  have  in  particular: 

{Cil3M  -h  A)(C'2M  -  A)(C3M  -  A)  X 
(A2  -  C4I3M)  =  0.  (59) 

It  is  seen  that  the  coefficients  Ci,  C2,  C3  and 
C4  scale,  respectively,  the  backward  acous¬ 
tic  speed,  the  shear-convection  speed,  the 
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entropy-convection  speed  and  the  forward 
acoustic  speed.  The  acoustic  wave  front, 
which  is  the  envelope  of  the  graph  of  the 
acoustic  plane  waves,  is  an  ellipse  with  equa¬ 
tion 

x-\{c,-c,)m\\  !  y 

i(C,  +  C4)  ) 

(60) 

In  order  to  make  the  elements  b  and  c  of  P 
equal  for  M  — *•  0,  we  must  satisfy 

^^  +  C,iM)=C4{M)f3  (61) 

for  M  — +  0,  or  simply 

Ci(0)  +  C2(0)  =  C4(0).  (62) 

This  means  that  both  the  backwaxd  acoustic 
speed  and  the  shear-convection  speed  must  be 
reduced  in  magnitude  with  respect  to  the  for¬ 
ward  acoustic  speed.  The  best  condition  num¬ 
ber  results  when  the  backward  acoustic  and 
shear  speeds  are  equal  in  magnitude,  i.  e. 

C:(0)  =  C2(0)  =  ic4(0).  (63) 

In  order  to  keep  the  normalization  that  the 
largest  characteristic  speed  equals  the  flow 
speed,  we  take  C4(0)  to  be  1.  which  makes 
both  Ci(0)  and  C2  equal  to  |.  Since  the  three 
coefficients  all  should  tend  to  1  for  the  larger 
subsonic  values  of  M,  there  is  no  reason  to  let 
C4  differ  from  1  for  any  Mach  number,  nor  Ci 
from  €2-  Our  final  choice  therefore  is 

Ci(M)  =  C2{M)  =  a{M),  (64) 


o{0)  =  i 

(65) 

a(l)  =  1, 

(66) 

C4(M)  =  1. 

(67) 

One  question  remains:  should  the  entropy- 
convection  speed  be  re-scaled,  or  may  it  re¬ 
main  equal  to  the  flow  speed?  There  is  some 
evidence  that  the  entropy  speed  better  equal 
the  shear  speed:  in  one  of  our  airfoil  calcula¬ 
tions  an  instabilty  encountered  at  the  far-field 
boundary  could  be  suppressed  by  lowering  tb 
entropy  speed  to  the  shear-speed  value.  A 
possible  explanation  is  that  for  the  precondi¬ 
tioned  system  of  equations  the  quantity  con- 
vected  besides  entropy  is  not  shear  strength 
but  actually  is  the  total  enthalpy  for  an  isen- 
tropic  process.  If  entropy  and  isentropic  to¬ 
tal  enthalpy  are  convected  at  different  speeds 
their  values  may  get  out  of  sync,  feeding  an 
instability.  If  this  explanation  is  correct,  an¬ 
other  way  of  preventing  the  instability  is  to 
add  a  multiple  of  the  entropy-convection  equa¬ 
tion  to  the  other  convection  equation,  so  that 
a  convection  equation  results  for  the  full  t( 
tal  enthalpy,  not  subject  to  a  speed  restric¬ 
tion.  However  this  may  be,  for  the  time  be¬ 
ing  we  shall  take  the  enthalpy-  and  entropy- 
convection  speeds  equal,  i.  e., 

C3(M)  =  a(M).  (68) 


Our  sub-optimal,  symmetric  preconditioner 
thus  becomes 


P(U)  = 


aMl 

P 

0 

0 

\ 

0 

a  (i  4- 1) 

0 

0 

0 

0 

/?. 

0 

0 

0 

0 

a 

/ 

(69) 
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For  our  numerical  experiments  we  chose 


a 

a 


a 


i  0<M<i 

2{i  +  3(m-1) 
1  „  2 

1, 


a  continuously  differentiable  function  which 
uses  a  cubic  to  switch  between  the  two  plateau 
values. 


4  Numerical  results 

We  experimented  with  both  preconditioners, 
the  cisymmetric  matrix  (52)  and  the  symmet¬ 
ric  matrix  (69);  the  results  can  be  summarized 
as  follows. 

1.  Both  new  asymmetric  and  symmetric 
preconditioners  made  it  possible  to  con¬ 
verge  with  a  conservative  upwind  code 
(structured  grid,  cell-average  update)  to 
steady  solutions  for  low-speed  flow  cases 
in  which  the  original  Van  Leer-Lee-Roe 
preconditioner  (8)  would  lead  to  instabil¬ 
ity  or  failure  to  converge.  The  symmet¬ 
ric  preconditioner  was  the  most  robust  of 
the  two,  allowing  the  use  of  lower  Mach 
numbers  and  finer  grids,  and  giving  faster 
convergence. 

2.  Even  the  better  preconditioner  of  the 
two  can  hardly  be  called  robust:  it  still 
needed  special  handling  in  order  to  avoid 


blow-up  for  an  impulsively  started  flow 
field,  in  addition  to  the  mandatory  down¬ 
ward  limiting  of  the  M- values  in  P  [2]. 

Surprisingly,  the  symmetric  precondi¬ 
tioner  greatly  improved  a  nonconserva¬ 
tive  flow  code,  namely,  the  unstructured 
cell-vertex  code  developed  by  Mesaros 

[13]  on  the  basis  of  fluctuation-splitting 
ideas.  The  equations  used  are  those 
found  after  preconditioning  the  Euler 
equations  with  the  Van  Leer- Lee- Roe 
preconditioning,  which  neatly  separates 
the  acoustic  equations  (a  pair  of  time- 
dependent  Cauchy-Riemann  equations) 
from  the  two  convection  equations  (one 
for  the  isentropic  total  enthalpy  and  one 
for  entropy).  The  convection  eqations 
are  updated  by  a  state-of-the-art  multi¬ 
dimensional  upwind  advection  scheme, 
the  acoustic  equations  are  by  a  cell- vertex 
distribution  scheme  similar  to  that  of  Ni 

[14] .  The  procedure  previously  failed  to 
converge  for  airfoil  flows  at  lower  inflow 
Mach-numbers,  with  the  strongest  inde¬ 
terminacy  occurring  in  the  leading-edge 
stagnation  region.  The  sub-optimal  sym¬ 
metric  preconditioner  decouples  the  equa¬ 
tions  just  as  the  optimal  one,  and  made  it 
possible  to  achieve  accurate  converged  re¬ 
sults  for  arbitrarily  low  M^o  on  a  fine  grid, 
without  any  further  special  measures. 

Figure  2  shows  Mach  contours  of  the 
steady  solutions  obtained  for  flow  over  a 
NACA  0012  at  zero  angle  of  attack,  for 
inflow  Mach  numbers  of  0.1  (top)  and 
0.01  (bottom).  The  grid  has  131  nodes  on 
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the  body,  strongly  clustered  at  the  lead¬ 
ing  edge.  The  two  solutions  are  practi¬ 
cally  indistinguishable  and,  as  shown  in 
Figure  1,  have  identical  convergence  his¬ 
tories,  except  for  a  scale  difference  of  a 
factor  100  between  the  residuals.  The  so¬ 
lutions  are  of  high  quality,  as  evidenced 
by  the  low  drag  coefficient  of  0.0002  ob¬ 
tained  in  both  cases. 

4.  When  changing  to  a  conservative  up¬ 
date  scheme  Mesaros’  code  lost  robust¬ 
ness  for  low-speed  flows,  showing  prob¬ 
lems  similar  to  those  experienced  with 
the  cell-average-based  finite- volume  code. 
In  particular,  the  code  would  not  toler¬ 
ate  impulsive-start  conditions;  the  non¬ 
conservative  version  of  the  code  had  to 
be  used  to  get  through  the  initial  tran¬ 
sients.  Alternatively,  increasing  the  ar¬ 
tificial  dissipation  in  the  acoustic  distri¬ 
bution  scheme  could  be  used  to  stabilize 
the  scheme  initially.  It  was  further  ob¬ 
served  that  the  lower  bound  on  M  in  the 
elements  of  P  had  to  be  raised  in  the  con¬ 
servative  code. 

5.  Another  conservative  cell-vertex  code 
[15],  though,  seemed  to  draw  no  benefit 
at  all  from  the  modified  preconditioner: 
it  performed  as  well  without  as  with  the 
modification.  This  code  was  developed  by 
D.  Darmofal  for  the  study  of  Euler  pre¬ 
conditioners,  in  particular,  the  effect  of  a 
preconditioner  on  the  eigenvector  struc¬ 
ture  of  the  equations.  The  scheme  imple¬ 
mented  was  Barth’s  [16],  which  uses  cell- 


vertex  values  to  define  Riemann  prob¬ 
lems  at  the  faces  of  co-volumes.  Dar¬ 
mofal  found  that  the  problems  of  insta¬ 
bility  or  nonconvergence  near  a  stagna¬ 
tion  point  axe  solely  caused  by  certain 
acoustic  eigenvectors  becoming  parallel 
for  M  — »  0,  and  that  bounding  M  away 
from  zero  in  the  (1,1)  element  of  P  solves 
these  problems.  So  far  the  following  com¬ 
binations  of  preconditioners  and  numeri¬ 
cal  flux  functions  have  been  tested: 

•  Turkel  preconditioning  with  scalar 
artificial  viscosity; 

•  Turkel  preconditioning  with  matrix 
viscosity  (upwind  differencing); 

•  Van  Leer  preconditioning  with 
scalar  axtificial  viscosity. 

Although  it  is  still  to  early  to  draw  a  firm 
conclusion,  it  seems  that  at  least  this  par¬ 
ticular  kind  of  cell- vertex  code  can  easily 
be  made  robust,  and  that  any  modifica¬ 
tion  of  the  preconditioning  should  be  me 
tivated  by  orthogonalizing  eigenvectors, 
rather  than  relaxing  the  condition  num¬ 
ber  of  eigenvalues. 

5  Conclusions 

In  the  present  study  we  deveoped  a  modified 
preconditioner  that  reduces  the  sensitivity  of 
the  preconditioned  flow  equations  to  the  flow 
angle  in  a  stagnation  region.  This  dramati¬ 
cally  improves  the  stability  and  convergence  of 
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at  least  one  flow  code,  a  nonconservative  cell- 
vertex  code  described  elsewhere  in  this  volume 
[13].  The  effect  on  certain  conservative  codes 
is  not  impressive,  while  yet  another  conser¬ 
vative  code,  also  described  elsewhere  in  this 
volume  [15],  does  not  seem  to  need  the  modi¬ 
fication  at  all.  All  preconditioners  need  some 
downward  limiting  of  the  value  of  the  Mach 
number;  its  beneficial  effect,  according  to  Dar- 
mofal  and  Schmid  [15],  is  mostly  to  prevent 
the  eigenvectors  of  certain  acoustic  waves  to 
become  parallel. 
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Figure  1:  Residual  histories  for  the  calcula¬ 
tions  of  the  steady  solutions  shown  in  Figure 
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Figure  2;  Mach  contours  for  flow  over  a  N  AC. A 
0012  at  zero  angle  of  attack;  for  a  descrip¬ 
tion  of  code  and  grid  see  the  main  text.  Top: 
Moo  =  bottom:  Moo  =  0.01. 
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Abstract 

An  optimization  formulation  is  described  for  multi¬ 
stage  schemes  based  on  the  hi-hi  high-frequency  content 
in  the  Fourier  footprint  of  the  preconditioned  spatial  op¬ 
erator.  These  coefficients,  when  used  in  conjunction  with 
semi-coarsened  multigrid  and  local  preconditioning  pro¬ 
vide  a  fast  and  robust  method  for  achieving  steady-state 
Euler  solutions.  Multigrid  speed-ups  of  3-4  times  are  ob¬ 
served  when  using  local  preconditioning  as  compared  to 
local  time-stepping.  The  extension  to  Navier-Stokes  op¬ 
erators  is  also  described. 

1  Introduction 

Explicit  marching  schemes  are  commonly  used  as 
multigrid  relaxation  schemes  when  solving  the  Euler  and 
Navier-Stokes  equations.  These  schemes  must  feature  ef¬ 
fective  high-frequency  damping  in  order  to  be  suited  for 
use  in  multigrid  marching.  Multi-stage  schemes  provide 
the  flexibility  to  achieve  the  desired  smoothing  proper¬ 
ties.  True  multigrid  convergence  rates  can  only  result  with 
multi-stage  relaxation  schemes  that  are  able  to  effectively 
damp  high-frequency  errors  for  all  flow  conditions. 

In  [1]  we  presented  an  optimization  scheme  to  obtain 
optimal  sequences  of  time-step  values  (i.e.,  a  multi-stage 
scheme)  for  discretizations  of  the  full  Euler  or  Navier- 
Stokes  spatial  operator.  ThougB  this  method  was  a  step 
forward  from  earlier  formulations,  which  based  the  de¬ 
sign  of  the  multi-stage  schemes  on  either  the  scalar  one¬ 
dimensional  [2,  3)  or  twodimensional  [4]  convection  equa¬ 
tion,  the  schemes  obtained  were  not  truly  independent  of 
flow  conditions,  such  as  Mach  number  and  flow  angle. 

In  [5]  Allmaras  pointed  out  that  with  semi-coarsening 
[6.  7],  the  high-frequency  domain  over  which  the  multi¬ 
stage  scheme  must  be  a  good  damper  of  errors  is  reduced 
to  “hi-hi”  combinations  (in  2-dimensions).  This  makes 
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it  easier  to  design  optimal  multi-stage  schemes  that  are 
largely  independent  of  flow  conditions.  In  [8]  we  have 
obteuned  such  multi-stage  schemes  and  have  shown  that 
these  multi-stage  schemes,  in  coi^unction  with  the  semi- 
coarsened  multigrid  algorithm  and  local  preconditioning 
[9],  provide  a  fast  and  robust  method  for  achieving  steady- 
state  Euler  solutions. 

In  this  paper  we  will  present  more  detailed  results, 
including  solutions  to  model  problems,  to  further  substan¬ 
tiate  this  claim.  Our  research  has  also  shown  that  Navier- 
Stokes  operators  require  a  different  optimization  formular 
tion  for  the  design  of  multi-stage  schemes.  This  issue  will 
also  be  addressed  in  this  paper. 


2  Semi-coarsening 

When  multi-dimensionaJ  convection  is  aligned  with 
one  of  the  grid  directions,  single-grid  relaxation  schemes 
cannot  damp  high-frequency  errors  propagating  in  the 
normal  direction  that  are  coupled  to  low  frequency  er¬ 
rors  in  the  convection  direction.  This  is  known  as  the 
single-grid  alignment  problem. 

Semi-coarsening  is  a  method  meant  to  resolve  this 
grid-alignment  problem  in  a  multigrid  context.  Mulder 
[6,  7]  has  developed  an  efficient  solver  for  the  steady  2-D 
Euler  equations  based  on  semi-coarsening.  The  method 
employs  semi-coarsening  in  two  directions  simultaneously 
(for  a  two-dimensional  problem). 

The  usual  restriction  and  prolongation  operators 
have  to  be  modified  to  handle  input  from  more  than  one 
grid.  If  one  grid  needs  data  from  two  finer  grids,  the 
two  sets  of  data  obtained  by  the  restriction  from  each 
finer  grid  are  averaged.  For  prolongation,  the  correction 
is  computed  with  respect  to  the  latest  fine-grid  solution 

As  pointed  out  earlier,  with  semi-coarsening,  the 
high-frequency  domain  over  which  the  multi-stage  scheme 
must  be  a  good  damper  of  errors  is  reduced  to  '‘hi-hf  com¬ 
binations,  making  it  easier  to  design  multi-stage  schemes 
that  are  largely  independent  of  flow  conditions. 
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Fourier  Footprint  of  hi-hi  domain. 

Modified  Roe  Scheme  with  M  =  0.1,^  =  0®,v  =  l. 
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Figure  1;  Fourier  footprint  of  the  first-order  upwind  ap¬ 
proximation  of  the  spatial  Euler  operator  (modified  Roe 
scheme)  with  the  preconditioner  of  Van  Leer  et  al.  [9],  for 
M  =  0.1,  and  flow  angle  (^  =  0*.  Footprint  corresponds 
to  hi-hi  domain  |^ri,  \0y  \  G  The  time-step  chosen 

corresponds  to  a  Courant-number  value  of  1. 

3  Local  Preconditioning 

Local  preconditioning  matrices  attempt  to  remove 
the  spread  among  characteristic  speeds  as  much  as  pos¬ 
sible.  The  matrix  derived  by  Van  Leer  et  al.  [9] 
achieves  what  can  be  shown  to  be  the  optimal  con¬ 
ditio  nnumber^f^^  characteristic  speeds,  namely, 

1  / ^1  -  min(M2,  where  Af  is  the  local  Mach  num¬ 

ber.  This  is  a  major  improvement  over  the  condi¬ 
tion  number  before  preconditioning,  which  equals  (A/  -h 
1)/  min(  Af,  \M  -  1|). 

The  effect  of  local  preconditioning  on  the  discretiza¬ 
tions  of  the  spatial  Euler  operator  is  a  strong  concentra¬ 
tion  of  the  pattern  of  eigenvalues  in  the  complex  plane. 
This  makes  it  possible  to  design  multi-stage  schemes  that 
systematically  damp  most  high-frequency  waves  admit¬ 
ted  by  a  particular  discrete  operator  [1].  The  resulting 
schemes  are  not  only  preferable  .as  solvers  in  a  multi-grid 
strategy,  they  are  also  superior  single-grid  schemes,  as  the 
preconditioning  itself  adrcady  accelerates  the  convergence 
to  a  steady  solution,  and  the  high-frequency  damping  pro¬ 
vides  robustness. 

The  local  preconditioning  matrix  described  in  [9]  was 
derived  based  on  the  partial  differential  equations  that 
make  up  the  Euler  equations.  On  analysis,  it  was  realized 
that  multigrid  damping  could  be  improved  by  modifying 
the  matrix  such  that  the  high-frequency  modes  of  all  the 
waves  overlap  more  completely.  This  modified  Euler  ma¬ 
trix  is  described  in  [10].  Figures  1  and  2  are  presented 
as  examples  of  the  hi-hi  frequency  content  in  the  Fourier 
footprints  obtained  with  this  modified  preconditioner. 
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Figure  2:  Fourier  footprint  of  the  k  =  1/3  upwind  ap¬ 
proximation  of  the  spatial  Euler  operator  with  the  pre¬ 
conditioner  of  Van  Leer  et  al.  [9],  for  Af  =  0.1,  and  flow 
angle  <f>  =  45®.  Footprint  corresponds  to  hi-hi  domain 
\0Ay  \0y\  €  (f ,  The  time-step  chosen  corresponds  to  a 
Courant-number  value  of  1. 

4  Optimization  procedure  for  Eu¬ 
ler  Operators 

The  procedure  for  optimizing  high-frequency  damp¬ 
ing  aims  at  minimizing  the  maximum  of  the  modulus  of 
the  scheme’s  amplification  factor  over  a  given  set  of  high- 
frequency  eigenvalues.  The  input  parameters  are  the  time- 
step  values  k  =  1, ..,  m,  of  an  m-stage  algorithm. 

When  updating  the  solution  of 

Ut  =  Res{U)  (1) 

from  time  level  t"  to  -f  At,  the  method  takes 

the  form 

=  f/",  (2) 

fj{k)  _  (7(0)4.  ,<:=  l,..,m,(3) 

fjn+l  _  ij{m)  (4) 

with  Af  =  Af^"*^.  According  to  linear  theory,  one  step 
with  the  full  scheme  multiplies  each  eigenvector  of  the 
operator  Res(U),  with  associated  eigenvalue  A,  by  a  factor 
of  the  form  „ 

TTl 

P{z)=  l-l-z-l-J^CkZ*,  (5) 

k-2 

where 

z  =  AAf  (6) 

generally  is  complex.  The  m  -  1  coefficients  c*  relate  to 
the  time-step  ratios  a*  =  Af^*V^^I  actual  time  step 
Af  is  the  mth  parameter. 
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Figure  3:  High-frequency  Fourier  footprint  of  the  precon¬ 
ditioned  first-order  upwind  Euler  operator  (modified  Roe 
scheme  with  preconditioner  of  Van  Leer  et  al.  plotted  on 
top  of  the  level  lines  of  the  amplification  factor  of  the  as¬ 
sociated  optimal  4-stage  scheme.  Flow  angle  0^,  M  =  0.9. 
Result  of  optimization  over  entire  high-frequency  domain 
minus  wedge  filter  (cf.  [1]). 

Given  a  discrete  Euler  operator,  the  optimization 
procedure  starts  out  by  computing,  for  a  fixed  combina¬ 
tion  of  M  and  (p  (=  flow  angle),  a  discrete  set  of  eigenval¬ 
ues  for  wave-number  pairs  in  the  high-frequency 

range,  i.e. 

|/?r|€(|,)r)  and  |/?y|  6  .  (7) 

Assuming  a  set  of  starting  values  for  the  m-stage 
scheme,  for  instance  Tai’s  values  [3],  the  value  of  |P(^)| 
is  computed  for  all  eigenvalues  previously  obtained, 
and  its  maximum  is  found.  This  is  our  functional 
Af,  o);  it  must  be  minimized  by  varying 
the  m  parameters. 

It  is  not  a  priori  clear  that  the  A<^*^  will  be  in¬ 
sensitive  to  values  of  M  auid  0.  If,  as  in  [1],  we  have 
to  consider  the  whole  high-frequency  domain,  i.e.  |/?r|  € 
(|,)r)  and/or  \0y\  6  two  problems  arise:  the 

alignment  problem  mentioned  earlier  and  a  singularity 
problem  for  iVf  — ►  1.  Figure  3  gives  evidence  of  both 
problems.  There  are  lo-hi  entropy/shear  modes  extending 
all  the  way  into  the  origin,  and  hi-lo  acoustic  error  modes 
at  some  distance  to  the  origin  that  vanish  as  \/l  - 
These  two  effects  go  away  if  we  restrict  ourselves  to  the 
high-frequency  range  (7)  appropriate  for  semi-coarsening. 

The  optimal  (in  the  Loo  sense)  m-stage  scheme  for 
a  given  discrete  Euler  operator  may  hence  be  obtained  as 
the  solution  to  the  following  minmax  problem: 

(Topt  =  min  (  max  l|P(-(/?rt  ;3y ,  i^),  5)|1  j  .  (8) 

(5,v)  ,r]  / 


Figure  4:  Optimal  four-stage  scheme  obtained  by  opti¬ 
mizing  over  hi-hi  frequency  footprint  of  the  modified  Roe 
scheme  with  the  preconditioner  of  Van  Leer  et  al.  [9]. 
M  =  0.1,  0  =  0.  aopt  =  0.0632,  i/  =  2.108. 

This  optimization  problem  is  solved  using  the 
method  of  simulated  annealing  in  conjunction  with  the 
downhill  simplex  algorithm  of  Nelder  and  Mead  [11,  12]. 

Simulated  annealing  has  proved  to  be  both  powerful 
and  robust  and  the  algorithm  does  not  require  frequent 
restarts  (unlike  Powell’s  method,  which  was  used  in  [1]). 

Figure  4  is  an  example  of  a  scheme  designed  with  the 
preconditioner  of  Van  Leer  et  al. 

4.1  Flow  angle  dependence 

Any  remaining  flow- angle  (0)  dependence  may  also 
be  removed  by  appropriate  definition  of  the  Courant  num¬ 
ber  u.  As  in  [1,  8],  for  the  preconditioned  Euler  equations, 
with  the  characteristic  speeds  equal  to  or  close  to  g,  we 
define  the  Courant  number  as 

i'  =  — IT’ 

/(Ax,  Ay,0) 

where  /  is  a  typical  cell-width  that  may  depend  on  the 
flow  direction.  For  rectangular  cells  we  find  that  defining 

/  =  l/(|cos(^l-l->^7ilsin<^|)  (10) 

takes  away  most  of  the  variation  of  u  with  flow  angle,  so 
that  a  single  value  may  be  recommended.  ([1]  has  this 
function  defined  incorrectly.  The  values  of  i/  in  the  tables 
in  this  reference  need  to  be  scaled  by  a  factor  of  \/2  when 
using  the  above  definition  of  /.  [13]  has  updated  tables). 

Figure  5  shows  the  variation  in  Courant  number  with 
length  scale  independent  of  flow  angle  for  a  typical  multi¬ 
stage  scheme.  As  can  be  seen  from  the  figure,  our  choice 
of  length  scale  provides  a  fairly  good  fit  to  the  variations. 
This  allows  us  to  recommend  a  single  value  of  v  for  each 
multi-stage  scheme. 
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Figure  5:  Variation  with  flow  angle  of  multi-stage  coeffi¬ 
cients  and  Courant  number  (based  on  a  fixed  length  scale 
independent  of  flow  angle)  for  a  first-order  4-stage  optimal 
scheme.  Optimization  over  hi-hi  frequency  domain  only. 


Figure  6:  Variation  of  multi-stage  coefficients  and  Courant 
number  with  Msich  number  for  a  first-order  4-8tage  opti¬ 
mal  scheme.  Optimization  over  hi-hi  frequency  domain 
only. 


4.2  Dependence  on  Mach  number 

The  Mach  number  dependence  of  the  multi-stage 
schemes  derived  in  the  previous  section  was  related  to  lo- 
hi  acoustic  error  modes  that  had  symbols  at  some  distance 
to  the  origin  that  vanish  as  Vl  -  A/^.  The  hi-hi  frequency 
error-modes  are  fairly  insensitive  to  Mach  number,  mak¬ 
ing  it  possible  to  design  optimal  multi-stage  schemes  with 
coefficients  that  are  fairly  insensitive  to  Mach  number  as 
well.  Figure  6  shows  the  variation  in  the  optimal  multi¬ 
stage  coefficients  with  Mach  number.  A  slight  change  in 
the  coefficients  is  observed  while  passing  through  the  sonic 
point,  but  the  change  in  values  is  small  enough  to  allow 
us  to  recommend  a  single  set  of  coefficients  for  each  mul¬ 
tistage  scheme  that  may  be  applied  over  the  entire  Mach 
number  and  flow-angle  range. 


5  Optimization  procedure  for 
Navier- Stokes  Operators 

For  Navier-Stokes  preconditionera  the  idea  put  forth 
in  [14]  is  to  make  the  sue  of  the- footprint  independent  of 
cell- Reynolds  number.  This  is  achieved  as  follows.  If  we 
write  the  2-D  discretized  Navier-Stokes  equations  as: 

U,  =  LeuU  -t-  (CU,)*  -I-  (DU,)v  -I-  (EUy)v,  (H) 

the  first  term  on  the  right-hand  side  is  the  discrete  Euler 
operator;  the  remaining  terms  are  the  viscous/conductive 
terms,  assumed  to  be  approximated  by  central  differenc¬ 
ing.  These  contribute  only  to  the  extent  of  the  footprint 
along  the  negative  real  axis,  which  is  inversely  propor¬ 
tional  to  the  cell-Reynolds  number.  The  proper  scaling 
required  to  make  the  size  of  the  footprint  independent  of 
cell-Reynolds  number  is  obtained  by  choosing  P/vs  as  fol¬ 


lows: 


p-1  _  p-1 
^  NS  —  ^ 


(12) 


For  higher-order  upwind  differencing,  the  same  sc2ding 
technique  for  the  highest  frequency  Fourier  footprints 
yields  a  similar  expression: 


P;;i  =  (l-«)Pii  +  ^c  +  iE.  (13) 


This  strategy  works  for  cell-Reynolds  numbers  that  are 
not  too  low,  e.g.  >0.1.  For  very  low  a  slight 

modification  is  needed  in  the  continuity  equation;  this  will 
be  discussed  elsewhere  (cf.  [15]). 

With  the  above  Navier-Stokes  preconditioning,  we 
are  able  to  design  a  family  of  multi-stage  schemes  that  are 
dependent  only  on  cell-Reynolds  number.  We  would  hope 
that  as  the  cell-Reynolds  number  decreases,  the  time  steps 
required  for  good  damping  increase,  as  the  high-frequency 
content  in  the  footprint  begins  to  align  itself  along  the 
negative  real  axis.  This,  unfortunately,  does  not  hap¬ 
pen  if  we  make  use  of  the  above  optimization  procedure 
for  a  discrete,  preconditioned  Navier-Stokes  operator  with 
prescribed  values  of  Mach  number,  flow  angle  and  cell- 
Reynolds  number.  The  formulation  described  above  at¬ 
tempts  to  maximize  the  damping  over  the  domain,  which 
results  in  extremely  (unnecessarily)  small  functional  (cr^pt) 
values  for  the  corresponding  multi-stage  schemes.  The 
stronger  damping  comes  at  a  price  -  smaller  Courant  num¬ 
bers  for  these  time-stepping  schemes.  This  means  taking 
smaller  time-steps,  thereby  increasing  the  number  of  com¬ 
putational  steps  required  to  attain  a  converged  solution. 
Since  the  cell  length-scales  used  in  computing  these  time- 
steps  tend  to  be  small  (a  common  feature  of  high- Reynolds 
number  computations),  the  larger  the  Courant  number  of 
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Figure  7:  Four-stage  scheme  obtained  by  optimizing  over  Figure  8:  Optimal  four-stage  scheme  obtained  by  opti- 
hi-hi  frequency  footprint  of  the  discrete,  preconditioned  mizing  over  hi-hi  frequency  footprint  of  the  discrete,  pre- 
Navier-Stokes  operator.  <Tmm  was  prescribed  as  0.1.  A/ =  conditioned  Navier-Stokes  operator.  M  =  0.1,  ^  =  0, 
0-1,  <^  =  0,  Recr  -  0.1,  =  1.  1/  =  7.891.  =  0.1,  AH  =  1.  (Topt  =  0.0034,  u  =  2.953. 


the  scheme,  the  better.  It  is  therefore  obvious  that  a  dif¬ 
ferent  approach  is  necessary  for  Navier-Stokes  operators. 

For  Navier-Stokes  operators,  we  have  redefined  the 
optimization  problem  as  follows:  For  a  given  spatial  oper¬ 
ator,  find  the  largest  Courant  number  with  which  the  cor¬ 
responding  multi-stage  scheme  has  a  prescribed  damping 
capability. 

This  may  be  written  as: 

Solve  (t(i/)  —  ^mirx  (14) 

where  is  defined  below: 

(r(i/)  =  min  I  max  ||^(^(/?r, /?y » 5)11  )  • 

The  solution  to  this  problem  is  fairly  simple. 
We  make  the  assumption  that  is  a  continuous, 

monotonically-increasing  function  within  our  range  of  in¬ 
terest.  (This  seems  to  be  a  valid  assumption  from  our 
experiments).  We  may  then  seeu'ch  for  a  root  to  equa¬ 
tion  (14)  using  an  algorithm  such  as  the  bisection  method 
(cf.  [11])  or  the  more  complex  algorithm  due  to  Brent 
[12]  which  combines  the  sureness  of  the  bisection  method 
with  the  speed  of  a  higher  order  method  when  appropri¬ 
ate.  Each  evaluation  of  cr(u)  requires  a  solution  to  the 
underlying  optimization  problem  (which  may  be  solved 
using  simulated  annealing  as  before).  Care  must  there¬ 
fore  be  taken  to  ensure  that  an  appropriate  set  of  starting 
values  of  a  are  used  for  each  evaluation  of  Another 

point  of  concern  is  the  choice  of  (Tmm-  Too  large  a  value 
of  (Tmtn  can  result  in  a  scheme  that  is  unstable  for  fre¬ 
quency  modes  other  than  the  “hi-hi”  combinations  being 
considered.  Figure  9  is  an  example  of  such  a  scheme. 


Figure  7  is  an  example  of  a  scheme  designed  using 
this  procedure.  In  contrast,  Figure  8  shows  a  scheme  de¬ 
signed  using  the  optimization  procedure  described  earlier. 
As  can  be  seen  from  the  figures,  the  two  schemes  have  dif¬ 
ferent  damping  properties.  By  prescribing  the  damping 
factor  (Tmin  to  be  0.1,  the  new  scheme’s  Courant  number 
increases  to  7.89  from  2.95  with  the  earlier  formulation. 

This  approach  is  not  well  suited  for  discretizations  of 
the  Euler  operator.  The  increase  in  Courant  number  is 
marginal  at  best.  Navier-Stokes  footprints  tend  to  align 
themselves  along  the  real  axis  in  the  complex  plane  with 
decreasing  cell-Reynolds  number  (as  the  behavior  becomes 
increasingly  parabolic).  Increasing  the  Courant  number 
causes  the  footprint  to  grow  primarily  in  this  direction 
for  these  operators.  Euler  footprints  grow  outward  in  all 
directions  at  the  same  rate,  making  it  more  diflBcult  to 
obtain  time-stepping  schemes  that  damp  all  error-modes. 

One  way  of  avoiding  multi-stage  coefficients  that  am¬ 
plify  some  error- modes,  as  in  Figure  9,  is  by  adding  in  a 
constraint  to  the  formulation.  An  appropriate  constraint 
would  be 

Pmat  =  max  |lP(z(^r,/?y,^),5)|l  <  1.  (15) 

1^3.1. I^»|€(0,t] 

This  constrained-optimization  problem  (equations 
(14)  and  (15))  can  be  reformulated  as  an  unconstrained- 
optimization  problem  by  making  use  of  a  penalty  function. 
We  can  write  ^{u)  as 

<t{v)  =  min(  max  \\P{z{0,c,  0y,v),S)\\ 

-b  7(/’ma*  -  1)').  (16) 

where  7  is  a  constant,  7  >  0.  The  larger  7  is,  the  more 
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Figure  9:  Four-stage  scheme  obtained  by  optimizing  over 
hi-hi  frequency  footprint  of  the  discrete,  preconditioned 
Navier-Stokes  operator.  <Tmin  was  prescribed  as  0.2.  M  = 
0.1,  =  0,  Re^r  =  QA.ATl  =  1.  i/  =  12.845.  It  is  obvious 

that  this  scheme  will  be  unstable  for  some  “hi-lo"’,  “lo-hi” 
and  “lo-lo”  frequency  combinations. 


Figure  10:  Four-stage  scheme  obtained  by  optimizing  over 
hi-hi  frequency  footprint  of  the  discrete,  preconditioned 
Navier-Stokes  operator  with  a  constraint  on  Pmax  intro¬ 
duced  via  a  penalty  function.  The  entire  footprint  is 
shown  here.  (Tmin  was  prescribed  as  0.2.  M  =  0.1,  ^  =  0, 
flcAx  =  0.1,  =  1.  i/  =  11.3679. 


strongly  is  the  constraint  enforced.  Taking  too  large  a 
value  of  7  however  can  affect  the  robustness  of  the  opti¬ 
mization  algorithm.  We  found  that  taking  taking  a  mod¬ 
est  value  of  7  initially  (around  100)  and  then  recomputing 
the  solution  with  this  trial  solution  and  7  around  1  x  10^ 
worked  well  in  enforcing  the  constraint.  Figure  10  is  an 
example  of  a  solution  obtained  in  this  manner.  We  have 
obtained  a  stable  scheme  with  only  a  modest  decrease  in 
Courant  number.  It  is  unclear  however  if  this  scheme 
will  perform  well  in  nonlinear  implementations.  Certain 
modes  are  likely  to  be  excited  if  the  time-step  taken  cor¬ 
responds  to  even  a  modest  increase  over  the  prescribed 
value.  One  option  is  to  modify  our  constraint  to  aillow  for 
some  variation  in  the  Courant  number  without  encounter¬ 
ing  amplification  factors  greater  thaui  one  for  some  wave- 
modes.  This  modified  constraint  could  be  written  as 

majc  f  max  l|^*(^i5)||) 


where  z  =  z{3i,3y,i/).  The  optimization  problem  may 
be  formulated  as  an  unconstrained  optimization  problem 
as  before.  Figure  11  is  an  example  of  the  solution  to  this 
problem  with  <5  =  0.2.  These  coefficients  are  likely  to  be 
more  reliable  in  a  multigrid  formulation  than  the  coeffi¬ 
cients  derived  earlier. 

Choosing  a  lower  value  for  (Tmim  say  around  0.12,  for 
the  case  considered  in  Figure  11,  would  obviate  the  need 
for  the  constraint  (eqn.  17)  and  its  associated  complica¬ 
tions. 

It  is  unclear  what  is  the  best  choice  for  (Tmin  in  terms 


of  obtaining  the  best  possible  multigrid  convergence  rates. 
This  could  likely  be  determined  from  multigrid  studies  in¬ 
volving  preconditioned  Navier-Stokes  operators  and  var¬ 
ious  test  conditions.  Multi-stage  schemes  with  =  0.2 
and  more  have  worked  well  in  our  multigrid  Euler  studies 
(though  we  did  not  attempt  to  prescribe  the  magnitude 
of  damping  here).  Similar  multigrid  studies  remain  to  be 
undertaken  with  the  preconditioned  Navier-Stokes  equa¬ 
tions. 


5.1  Multi-stage  coefficients  with  prescribed  damp¬ 
ing  for  Navier-Stokes  operators 

When  following  this  approach,  it  is  unlikely  that  we 
would  be  able  to  recommend  a  single  set  of  coefficients 
independent  of  Re^x  and  A  likely  compromise  would 
be  to  generate  curve-fits  to  account  for  the  variations  in 
coefficients  due  to  each  of  these  parameters. 

As  a  preliminary  study  we  have  obtained  the  coef¬ 
ficients  for  a  four-stage  time  stepping  scheme  based  on 
the  Fourier  footprint  of  the  discrete,  preconditioned  first- 
order  Navier-Stokes  operator,  discretized  using  the  mod¬ 
ified  Roe  scheme  for  inviscid  terms  and  central  differenc¬ 
ing  for  the  viscous  terms.  The  variation  in  the  coefficients 
with  cell-Reynolds  number  was  studied.  Figure  12  shows 
these  variations,  normalized  to  the  values  of  the  coeffi¬ 
cients  for  the  Euler  operator.  The  figure  indicates  that 
the  optimal  Courant  number  increases  sharply  as  the  cell- 
Reynolds  number  decreases,  while  the  multi-stage  coeffi¬ 
cients  decrease.  It  would  seem  that  curve-fits  in  terms  of 
log(/ZeAr)  are  feasible. 
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Figure  11:  Four-stage  scheme  obtained  by  optimizing  over 
hi-hi  frequency  footprint  of  the  discrete,  preconditioned 
Navier-Stokes  operator  with  a  constraint  on  (eqn. 

17)  introduced  via  a  penalty  function.  The  entire  footprint 
is  shown  here.  <Tmin  was  prescribed  as  0.2.  M  =  0.1, 
=  =  0.1,  A1l=  I,  6-  0.2.  u  =  9.5361. 


Figure  12:  Variation  with  cell-Reynolds  number  of  mul¬ 
tistage  coefficients  and  Courant  numver  for  a  first  order 
4-stage  scheme  with  prescribed  damping.  Optimization 
based  hi-hi  frequency  footprint  of  the  discrete,  precondi¬ 
tioned  Navier-Stokes  operator  with  a  constraint  on  Pmax 
introduced  via  a  penalty  function.  (Tmin  was  prescribed  as 
0,15.  M  =  0.1,  0  =  0,  =  1,  5  =  0.05. 


5.2  Optimal  multi-stage  schemes  for  the  precon¬ 
ditioned  Euler  equations 

We  have  computed  optimal  multi-stage  schemes 
based  on  the  modified  Roe  discrete  operator  for  the  pre¬ 
conditioned  Euler  equations  with  the  preconditioner  of 
Van  Leer  et  al.  These  coefficients  are  presented  as  ta¬ 
bles  in  [13,  8],  and  have  been  used  to  obtain  the  results 
in  this  paper.  For  completeness,  a  condensed  version  of 
these  tables  is  presented  here;  see  Tables  1-4, 

These  schemes  are  not  only  preferable  as  solvers 
in  a  multigrid  strategy,  but  are  also  superior  single-grid 
schemes,  as  the  preconditioning  itself  already  acceler¬ 
ates  the  convergence  to  a  steady  solution  and  the  high- 
frequency  damping  provides  robustness. 

6  Multigrid  Euler  Solutions 

The  preconditioner  of  Van  Leer  et  al.  has  been  shown 
to  produce  dramatic  speed-ups  for  steaidy  solutions  to  the 
problem  of  '‘flow  past  a  semi-circular  bump  (t/c  =  0.042) 
in  a  channel”  [8].  The  modified  Roe  scheme  that  is  used 
here  also  provides  more  accurate  solutions  for  cases  with 
low  freestream  Mach  number.  Results  of  these  runs  are 
shown  in  Tables  5  and  6  in  terms  of  equivalent  work  units. 
A  work  unit  is  the  amount  of  work  required  to  compute  a 
single-stage  update  on  the  fine  grid.  These  are  the  “best” 
single  and  multigrid  results  obtained  from  multiple  runs 
using  2-6  (3-6  for  second  order)  stages  and  1-5  (1-4  for  sec¬ 
ond  order)  grid  levels.  The  second-order  solutions  made 
use  of  Van  Albada’s  limiter  and  defect-correction  multi- 
grid  cycles  for  the  multigrid  cases.  Figures  14,  15  and  16 


are  examples  of  solutions  obtadned  with  this  method. 

Tables  5  amd  6  indicate  that  local  and  matrix  time¬ 
stepping  require  similar  amounts  of  work  in  the  single-grid 
subsonic  and  transonic  runs.  This  is  probably  an  artifact 
of  the  problem.  (This  behavior  is  not  observed  when  solv¬ 
ing  for  the  flow  in  a  semi-infinite  channel).  The  reflective 
wall-boundary  conditions  in  the  channel  provide  minimal 
attenuation  of  acoustic  error-modes  and  this  effect  dom¬ 
inates  single-grid  convergence.  This  is  not  a  problem  for 
the  supersonic  cases,  which  have  a  largely  convective  na¬ 
ture.  The  speed-up  of  multigrid  with  local  precondition¬ 
ing  over  multigrid  with  local  time-stepping  is  a  lot  more 
dramatic;  the  former  is  3-4  times  faster  in  all  cases.  Con¬ 
vergence  in  the  subsonic  and  transonic  cases  is  also  an  or¬ 
der  of  magnitude  faster  than  in  the  corresponding  single¬ 
grid  cases.  Local  preconditioning  performs  admirably  on 
a  single-grid  for  the  supersonic  case,  and  there  is  not  much 
improvement  possible  without  modifying  multigrid  to  bet¬ 
ter  handle  convection-dominated  flow. 

A  semi-coarsened  multigrid  cycle  is  more  expensive 
than  a  regular  multigrid  cycle  (a  factor  of  three  on  aver¬ 
age  for  the  cases  considered  -  this  factor  depends  on  the 
number  of  grid  levels  and  number  of  stages  chosen).  Each 
semi- coarsened  multigrid  cycle  is  able  to  reduce  the  resid¬ 
ual  norm  by  a  larger  factor  than  a  corresponding  regular 
multigrid  cycle.  However,  some  of  the  extra  overhead  that 
is  inherent  in  the  method  does  show  up  in  the  work  re¬ 
quired  to  obtain  the  same  solution  with  the  same  number 
of  grid  levels  and  stages  in  the  time-stepping  scheme.  Ta¬ 
bles  5  and  7  refect  the  difference  in  work  required.  It 
should  be  pointed  out  though  that  the  semi-coarsened 


7 


t/c  =  0.042,  modified  Roc  scheme  with  «  =  0. 


Figure  14:  Mach  number  contours.  Subsonic  flow  past 
a  bump  in  a  channel.  Modified  Roe  scheme  with  the 
preconditioner  of  Van  Leer  et  al.  «  =  0  variable  ex¬ 
trapolation.  128x64  grid.  Solution  obtained  with  semi- 
coarsened  multigrid  incorporating  an  optimal  four-stage 
scheme.  M  =  0.05,  t/c  =  0.042. 


Figure  13:  Residual  history  plots  for  the  transonic  test 
case  considered  above  (Af  =  0.85).  64x32  grid.  All  runs 
were  made  with  four-stage  time-stepping  schemes  and  all 
multigrid  runs  made  with  four  levels  (8x4  as  coarsest  grid). 

multigrid  method  makes  up  for  the  extra  work  required  in 
the  form  of  increased  robustness.  Residual  history  plots 
indicate  that  the  convergence  rate  with  regular  multigrid 
can  be  erratic,  especially  in  the  transonic  regime,  whereas 
it  is  essentially  constant  with  semi-coarsening  (see  Fig¬ 
ure  13).  Also,  in  our  calculations,  regular  multigrid  ex¬ 
hibited  extremely  poor  convergence  rates  (especially  with 
local-timestepping)  when  attempting  to  obtain  a  second- 
order  spatially  accurate  solution  to  the  transonic  case 
{\f  =  0.85). 

6.1  Improving  robustness 

The  preconditioner  of  Van  Leer  et  al.  is  known  to  be 
lacking  in  robustness  about  sonic  and  stagnation  points. 
The  stagnation-point  singularity  appears  to  be  quite  a  dif¬ 
ficult  problem  to  correct  and  there  has  been  substantial 
research  in  this  area  (cf.  [16,  17]). 


f/c  =  0.042,  modified  Roe  scheme  with  «  =  0. 

Imn  0.618 
2  0.680 
4'0.760 
6  0.840 
8  0.920 
10  1.000 
12  1.080 
14  1.160 
16  1.240 
mx  1.293 

""T . . . 

3.0 

Figure  15:  Mach  number  contours.  Transonic  flow  past  a 
bump  in  a  channel.  Modified  Roe  scheme  with  the  pre- 
conditioner  of  Van  Leer  et  al.  k  —  0  variable  extrapo¬ 
lation  with  Van  Albada’s  limiter.  128x64  grid.  Solution 
obtained  with  semi- coarsened  multigrid  incorporating  an 
optimal  four-stage  scheme.  M  =  0.85,  t/c  =  0.042. 


t/c  =  0  042,  modified  Roe  scheme  with  x  =  0. 


mn  I  039 
2  1  081 
4  1  164 
6  1  247 
8  1  330 
10  1  413 
12  1  495 
14  1  578 
mx  1  660 


We  opted  to  fit  a  parabolic  curve  to  /?  such  that  the 
fit  matched  the  curve  \/l  -  for  a  given  value  of  M  in 
both  value  and  slope.  The  point  of  intersection  with  the 
other  branch,  \/M^  -  1,  was  also  required  to  match  in 
slope.  Since  only  three  conditions  can  be  specified  for  this 
curve-fit,  we  opted  not  to  specify  the  point  of  intersection 
with  the  curve  -  1  (cf.  [13]). 

With  this  approach  and  in  conjunction  with  explicit 


Figure  16:  Mach  number  contours.  Supersonic  flow  past 
a  bump  in  a  channel.  Modified  Roe  scheme  with  the  pre¬ 
conditioner  of  Van  Leer  et  al.  k  =  0  variable  extrapo¬ 
lation  with  Van  Albada’s  limiter.  128x64  grid.  Solution 
obtained  with  semi-coarsened  multigrid  incorporating  an 
optimal  four-stage  scheme.  M  =  1.4,  t/c  =  0.042. 
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Number  of  Stages 


2 

3 

4 

5 

6 

B 

0.3333 

0.1467 

0.08125 

1 

0.3979 

0.2033 

0.1240 

0.08516 

0^3 

1 

0.4226 

0.2343 

0.1521 

m 

1 

0.4381 

0.2566 

m 

1 

0.4525 

0^6 

1 

Bl 

1.0000 

1.5252 

2.1058 

2.6824 

3.0827 

^opt 

0.3333 

0.1418 

0.06328 

0.03024 

0.01627 

Table  1:  Optimal  multi-stage  coefficients  for  first-order 
scheme.  Optimization  based  on  hi-hi  high  frequency  do¬ 
main. 


Number  of  Stages 


:Hl 

2 

3 

4 

5 

6 

0.6826 

0.2695 

0.07727 

aj 

1 

1 

0.2428 

Qf4 

1 

0.5887 

0.3650 

OCB 

1 

0.5729 

Ole 

1 

0.7167 

1.9924 

2.4407 

O'opt 

0.6677 

0.4454 

0.3171 

0.2336 

0.1686 

Table  4:  Optimal  multi-stage  coefficients  for  /c  =  1/3 
scheme.  Optimization  based  on  hi-hi  high  frequency  do¬ 
main. 


Number  of  Stages 


2 

3 

4 

5 

6 

Qfl 

0.5713 

0.2239 

0.1299 

0.08699 

0.06134 

1 

0.5653 

0.2940 

0.1892 

0.1322 

1 

0.5604 

0.3263 

0.2201 

Q4 

1 

0.5558 

0.3425 

05 

1 

0.5531 

06 

1 

1/ 

0.6305 

1.0458 

1.4008 

1.7471 

2.0701 

O'opt 

0.6475 

0.4279 

0.2927 

0.2047 

0.1464 

t/e  =  0.042,  modified  Roe  scheme  with  =  0  and  expUct  RS. 


Fmin  0.63S 
Fst  0.640 
Fine  0.030 
Fmax  1.404 


Figure  17:  Mach  number  contours.  Flow  past  a  non¬ 
smooth  bump  in  a  channel.  Af  =  1.0,  ije  =  0.042, 
128  X  64  grid,  k  =  0  variable  extrapolation  with  modified 
Roe  scheme  and  Van  Albada’s  limiter.  Explicit  residual- 
smoothing  was  used  to  obtain  this  solution. 


Table  2:  Optimal  multi-stage  coefficients  for  k  =  0 
scheme.  Optimization  based  on  hi-hi  high  frequency  do¬ 
main. 


Number  of  Stages 


2 

3 

4 

5 

6 

Oil 

0.4450 

0.1780 

0.09900 

0.06431 

0.04540 

Q2 

1 

0.4774 

0.2434 

0.1509 

0.1044 

<^3 

1 

0.4913 

0.2783 

0.1846 

Ol4 

1 

0.5004 

0.3030 

05 

1 

0.5116 

05 

1 

1/ 

0.4386 

0.7439 

1.0139 

1.2608 

1.4613 

O’ opt 

0.6154 

0.3916 

0.2526 

0.1654 

0.1145 

Table  3:  Optimal  multi-stage  coefficients  for  /c  =  -1 
scheme.  Optimization  based  on  hi-hi  high  frequency  do- 


residual  smoothing,  we  were  able  to  obtain  a  solution  to 
the  channel  problem  with  a  freestream  Mach  number  .V/  = 
1.  Figure  17  shows  this  solution. 

7  Multigrid  for  Navier-Stokes  Op¬ 
erators 

Just  as  for  Euler  discretizations,  the  multi-stage 
schemes  described  above  for  discrete  Navier-Stokes  oper¬ 
ators  should  be  an  ideal  choice  as  the  basic  relaxation 
scheme  in  a  semi-coarsened  multigrid  strategy.  Multignd 
applications  of  the  optimized  multi-stage  schemes  will  be 
presented  in  [15]. 

8  Conclusions 

We  have  demonstrated  that  the  combination  of  local 
preconditioning  and  multi-stage  time-stepping  can  pro¬ 
duce  relaxation  schemes  that  boast  guaranteed,  strong 


main. 
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M  =  0.35 

M  =  0.85 

M  =  1.4 

Local  TS 

Matrix  TS 

Single 

2940 

2520 

4140 

II^QIII 

1550 

355 

Multigrid 

1123 

268 

911 

326 

651 

234 

Table  5:  Comparison  of  first-order  convergence  rates  for  flow  past  a  semi-circular  bump  in  a  channel,  64x32  grid. 
Work  units  required  to  reduce  the  norm  of  the  residual  by  five  orders  of  magnitude.  (Local  TS  =  local  time-stepping, 
Matrix  TS  =  local  preconditioning).  Semi-coarsened  multigrid. 


II 

o 

M  =  0.85 

M  =  1.4 

Local  TS 

Matrix  TS 

Local  TS 

Matrix  TS 

Local  TS 

Matrix  TS 

Single 

3685 

2380 

5305 

1344 

414 

Multigrid 

582 

191 

515 

167 

722 

309 

Table  6;  Comparison  of  second-order  («  =  0)  convergence  rates  for  flow  past  a  semi-circular  bump  in  a  channel,  64x32 
grid.  Work  units  required  to  reduce  HT^Hi  (cf.  [7])  to  10”^  for  M  =  0.35  and  M  =  1.4  and  to  5  x  10  ^  for  M  —  0.85. 
Nested  iteration  with  5  defect-correction  sweeps  on  each  coarse  grid  level  was  used  initially  to  improve  robustness  for 
multigrid  solutions.  Semi-coarsened  multigrid. 


M  =  0.35 

M  =  0.85 

M  =  1.4 

Local  TS 

Matrix  TS 

Local  TS 

Matrix  TS 

Local  TS 

Matrix  TS 

Single 

288 

Multigrid 

165 

174 

161 

Table  7:  Comparison  of  first-order  convergence  rates  for  flow  past  a  semi-circular  bump  in  a  channel,  64x32  grid 
Work  units  required  to  reduce  the  norm  of  the  residual  by  five  orders  of  magnitude.  Regular  multigrid. 
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high  frequency  damping  for  the  entire  range  of  flow  angles, 
Mach  numbers  and  cell-Reynolds  numbers.  Such  schemes 
are  ideally  suited  for  use  u  single-grid  relaxation  schemes 
in  a  multigrid  relaxation  framework,  particularly  if  semi¬ 
coarsening  is  used. 

Some  numerical  results  are  presented  to  support  this 
claim;  more  multigrid  experimenting  is  needed  to  deter¬ 
mine  the  best  balance  between  parameters  such  as  damp¬ 
ing  rate  and  time-step  value. 
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